Mountains, climate and niche heterogeneity explain global patterns of fern diversity
Editor: Gerald Schneeweiss
Abstract
Aim
It is well known that the distribution of species diversity is spatially heterogeneous, but understanding the factors contributing to this heterogeneity and to the formation of biodiversity hotspots remains a challenge. Here, we seek to improve our understanding of how historical, ecological and evolutionary processes contribute to current patterns of global fern diversity.
Location
Worldwide.
Taxon
Ferns.
Methods
To evaluate the drivers of global fern diversity, we integrate over 800,000 georeferenced species occurrence records of nearly 8000 species, a time‐calibrated phylogeny and seven climate and environmental layers. We use these data to summarize diversity and evolutionary patterns at a resolution of 100 × 100 km, and identify hotspots of fern species richness and endemism. We compare these hotspots to neighbouring non‐hotspot regions to provide insight into the factors controlling global patterns of fern diversity.
Results
Tropical and subtropical mountains harbour a disproportionate amount of species relative to the land area they occupy; 58% of global species richness occur in eight principally montane hotspots together comprising just 7% of Earth’s land area. We identify hotspots of fern species richness and endemism that are universally characterized by disproportionately high ecological variation. We demonstrate that total fern species richness scales linearly with available climate space at regional and global scales.
Main Conclusions
Areas of high environmental heterogeneity harbour a disproportionate amount of fern species, and global patterns of extant fern diversity reflect the distribution of these areas, especially in mountains at lower latitudes. Persistence of ancient lineages in areas with long‐term climatic stability helps explain exceptional endemism in regions such as Malesia.
INTRODUCTION
Tropical mountains harbour a disproportionate amount of Earth's diversity compared to adjacent lowlands and temperate regions, and are frequently recognized as biodiversity hotspots (Antonelli et al., 2018; Hoorn et al., 2013; Körner, 2000; Rahbek, Borregaard, Antonelli, et al., 2019). Although the mechanisms driving mountain biodiversity remain poorly understood, recent advancements indicate the interplay of several factors including, the availability of diverse niche space, habitat generation and ecological stability through time may contribute to these patterns (Perrigo et al., 2020; Rahbek, Borregaard, Colwell, et al., 2019). Mountains and other topographically complex regions in the tropics harbour an exceptional density of climatic niche space, compared to adjacent lowlands and temperate mountains (Rahbek, Borregaard, Antonelli, et al., 2019). Rapid mountain uplift also generates new habitats into which species can disperse to and diversify (Antonelli et al., 2018; Hoorn et al., 2013; Xing & Ree, 2017). Additionally, mountains and other climatically isolated regions can act as refugia for lineages by harbouring stable niches for lowland or cold‐adapted species to inhabit during climatic fluctuations over geologic time (Fjeldsa et al., 1999; Schönswetter et al., 2005; Tang et al., 2018). Teasing apart the contributions of these processes and understanding their interactions are critical for developing a cohesive understanding of the distribution of the world’s biota.
Biodiversity hotspots are maintained through a variety of mechanisms, including exceptionally high speciation, low extinction, immigration or a combination of these factors (Harrison & Noss, 2017; Jablonski et al., 2006; McKenna & Farrell, 2006; Rahbek, Borregaard, Antonelli, et al., 2019; Wiens & Donoghue, 2004). Along this spectrum, regions have been labelled as either ‘cradles’ or ‘museums’ of biodiversity in reference to places where speciation is high and extinction is low (Haffer, 1969; Richardson et al., 2001), or regions where extinction is low and immigration may be high, respectively (Fischer, 1960; Gaston & Blackburn, 1996; Stebbins, 1974). Prominent examples of cradles include montane regions that have experienced recent rapid uplift such as the tropical Andes and the Hengduan Mountains (Hoorn et al., 2010; Lagomarsino et al., 2016; Testo et al., 2019; Xing & Ree, 2017). Steep tropical or subtropical mountains like these are climatically stratified and seasonal variation in temperature is low compared to temperate mountains (Janzen, 1967). The interaction between high climatic heterogeneity across strata and strong seasonal stability within strata is hypothesized to lead to a decrease in species range size and niche breadth, limiting dispersal and gene flow, and promoting speciation (Janzen, 1967; Polato et al., 2018). Older mountains and other geologically or climatically stable regions, in contrast, have been demonstrated to act as refugia (Cronk, 1997; Fjeldsa et al., 1999; Tang et al., 2018). Broad macroevolutionary studies incorporating phylogenetic, ecological, climatic, topographic and georeferenced data are needed to untangle these processes to better understand the distribution of life on Earth, but until recently these studies have remained sparse due to limitations of data and the difficulty of modelling these processes (Mishler et al., 2014; Moreau & Bell, 2013).
We attempt to understand the factors underlying the exceptional richness of tropical montane biota by looking at ferns, a species‐rich group of vascular plants whose diversity peaks in tropical montane ecosystems (Kreft et al., 2010; Moran, 1995, 2008). Ferns are a compelling study group for addressing these questions because they are diverse (> 11,000 species; PPG, 2016), have a densely sampled and well‐resolved phylogeny (Testo & Sundue, 2016) and exhibit marked variation in their richness at regional and global scales, which is thought to be driven principally by eco‐climatic factors such as precipitation regimes and edaphic specialization (Khine et al., 2019; Kreft et al., 2010; Lehtonen et al., 2017; Salazar et al., 2015; Sundue et al., 2015; Weigand et al., 2020), rather than biotic interactions, as they produce wind‐dispersed spores (Kluge & Kessler, 2007). Although the origin of the fern lineage dates to the Devonian, much of the extant diversity accrued relatively recently (Schneider et al., 2004). This mix of persistent depauperon lineages with recent radiations provides an exceptional opportunity to study the accumulation of diversity and underlying ecological factors across space and time. Despite the group’s prominence, there have been few attempts to quantify fern diversity at a global scale (Christ, 1910; Kessler et al., 2011; Lehtonen et al., 2017; Weigand et al., 2020). To establish a global perspective, we combine a time‐calibrated molecular phylogeny and more than 800,000 cleaned GBIF occurrence records of almost 8000 species to identify hotspots of fern species richness, endemism and speciation. We integrate this framework with climatic, topographic and edaphic data to understand which factors underlie global and regional variation in fern diversity. We ask the following questions: Where are fern diversity hotspots? Do diversity hotspots occupy a disproportionate climatic and ecological niche space compared to adjacent regions? Do ferns speciate more rapidly in areas with high ecological space compared to more climatically stable regions? Do montane components of tropical hotspots harbour exceptional diversity compared to lowland regions in the same hotspot?
1 MATERIALS AND METHODS
1.1 Speciation rate estimation
Lineage speciation rates were estimated using a time‐calibrated phylogeny of 3973 fern taxa (Testo & Sundue, 2016), implementing a branch‐specific birth–death model in RevBayes (Höhna et al., 2016). We accounted for incomplete taxon sampling by assigning family‐level sampling proportions, using species richness estimates from the Pteridophyte Phylogeny Group (PPG, 2016). We specified a prior on the number of rate categories at k = 6 and assigned lognormal distributions for speciation and extinction rates with means equal to 4.64 ((number of taxa in phylogeny/root age)/2) and standard deviations set to 0.5. We allowed for rate shifts to occur across the tree and set a uniform distribution on the number of shifts scaled between 0 and 10. This model was sampled with an MCMC on the University of Florida’s HiPerGator 2.0 high‐performance computing cluster for 24,000 generations, with 4000 generations specified as burn‐in. Convergence and mixing were assessed using Tracer 1.7 (Rambaut et al., 2018) to evaluate run convergence and to ensure that all parameters had an effective sample size of > 200. Tip speciation rates were extracted and used for subsequent analyses.
1.2 Spatial analyses
Species occurrence records were downloaded from GBIF with queries filtered to include only georeferenced fern herbarium specimens, resulting in a total of 1,686,374 records (http://GBIF.org“ GBIF.org (22 April 2020) GBIF Occurrence Download https://doi.org/10.15468/dl.amnnxg”.). To improve the quality of these records, we carried out a multi‐step occurrence cleaning process. First, we used the R package ‘CoordinateCleaner’ (Zizka et al., 2019) to flag and remove records that satisfied any of the following criteria: 0,0 coordinates, coordinates in an ocean, coordinates within 5 km of country centroids and capitals, records within 1 km of biodiversity institutions, and records with reversed latitude and longitude values; this removed 809,228 records. After removing these geographical errors, we resolved the taxonomy by passing all remaining records through a taxonomic thesaurus developed from the Checklist of Ferns and Lycophytes of the World (https://worldplants.webarchiv.kit.edu/ferns/), with additional inputs from Tropicos (www.tropicos.org) and the Pteridophyte Collections Consortium data portal (www.pteridoportal.org). The taxonomic thesaurus included a total of 59,448 published names (10,657 accepted names). All records were matched against this list and records with names that were classified as synonyms were updated to the corresponding accepted name; this resulted in name changes for 197,875 records. A total of 2386 records could not be matched to an accepted name and were excluded. Next, we used the R package ‘rbokeh’ (Hafen & Continuum Analytics, 2016) to manually generate interactive maps of species occurrences (visualized at the genus level) to identify further erroneous records, species outside of their native range and flag suspicious records for additional inspection; an additional 29,299 records were manually removed, resulting in a total of 845,878 records representing 7865 total species.
We summarized regional fern diversity by associating cleaned occurrence records with 100 × 100 km grid cells using the ‘join attributes by location’ tool in QGIS 3.10.1 and trimming the resulting data to retain distinct species per grid cell. For each grid cell, we calculated species richness and weighted endemism (Williams & Humphries, 1994), which represents the mean of the inverse range sizes of the taxa represented in each grid cell weighted by species richness. Weighted endemism aims to avoid issues with the use of arbitrary range sizes used to define an endemic taxon (Slatyer et al., 2007). All maps were depicted with an equal‐area (Mollweide) projection.
Hotspots of species richness and endemism were identified using the Gi* local statistic (Getis & Ord, 1996) implemented in QGIS with the Hotspot Analysis plugin (Oxoli et al., 2017). This approach allows for the detection of spatial autocorrelation of a given variable(s) (in this case, species richness and endemism, which were calculated for each grid cell), thereby allowing clusters of user‐defined areas (in this case, 100 × 100 km grid cells) of high values to be identified. We allowed for distance band optimization with a distance step of 10 km and identified global fern diversity hotspots as regions of adjacent grid cells with Z scores ≥ 1.65 for either species richness and endemism or both. Once hotspot cells were identified, we also designated nearby ‘non‐hotspot’ areas in three broad regions for comparative purposes. These non‐hotspots included all grid cells that were not identified as hotspots of either richness or endemism in the following areas: (1) all of continental America from 23ºN to 35ºS (tropical and southern subtropical America) and (2) all of Africa south of 10ºN (sub‐Saharan Africa), and Asia/Malesia from 70–155ºE and 30ºN–10ºS (eastern and southern Asia). Species richness and endemism estimates were obtained by extracting occurrence data within the cells comprising each hotspot; species were considered endemic to a hotspot if no occurrence records for the species were recovered outside of that hotspot. To generate estimates of species richness per elevational band, we extracted elevation values from the GMTED DEM layer (see ‘Section 2.3’) for each occurrence point and then sorted records by that value.
Mean speciation rates of lineage assemblages in each grid cell (average grid cell speciation rates) were calculated based on the 3218 species for which rates were estimated (these were species represented both in our occurrence data and the dated phylogeny of Testo & Sundue, 2016). A total of 604,748 occurrence points for species with speciation rate data were obtained from our total dataset and a single representative per grid cell was retained to calculate average grid cell speciation rates. Using these occurrences and time‐calibrated phylogeny, we also calculated spatial phylogenetic indices of phylogenetic endemism (PE) and phylogenetic diversity (PD) for each grid cell. PD is the sum of the total branch lengths that connect sets of taxa to a most recent common ancestor (Faith, 1992). PE is the range‐size‐weighted equivalent of PD (Laffan et al., 2016; Rosauer et al., 2009). These indices were calculated for use in downstream analyses, which are described below.
We identified centres of neo‐ and palaeo‐endemism with the CANAPE analysis (Mishler et al., 2014), which identifies areas with different endemism signals through a two‐step process of evaluating relative PE (RPE) values. RPE is the ratio comparing the observed PE from the original tree to PE where all branches of the tree are of equal length (PE_alt). In the first step, we used Biodiverse 3.0.0 (Laffan et al., 2010) to calculate PE, PE_alt and RPE values for all grid cells in our dataset and identify grid cells with significantly high values for either PE or PE_alt; significance (two‐tailed test, α = 0.05) was assessed by comparison to a null distribution of RPE values generated with 999 iterations of the ‘rand_structured’ function. In the second step, grid cells with significantly high values of either PE or PE_alt (or both) were then divided into four endemism classes based on their RPE scores. Grid cells with significantly high RPE are considered to be centres of palaeo‐endemism, cells with significantly low RPE are considered to be centres of neo‐endemism, cells with non‐significant RPE but significantly high PE and PE_alt are considered mixed endemism sites; in cases where both PE and PE_alt are significant at the α = 0.001 level, the cell is considered super‐endemic.
1.3 Environmental and climatic data
To compare environmental niche space of diversity hotspots and non‐hotspots, we obtained climate, soil and topographic data from CHELSA (Karger et al., 2017), SoilGrids (Hengl et al., 2017) and the USGS GMTED2010 (Danielson & Gesch, 2011) databases. The variables used were as follows: annual mean temperature, minimum temperature of the coldest month, annual precipitation, precipitation seasonality, soil pH, soil carbon content and mean elevation. Climate and elevation layers used were at a 30‐arc‐second resolution; soil variables were at a 15 cm depth and a 250 m resolution. Values for each variable were extracted at specimen occurrence points using the ‘sample raster values’ function in QGIS; values were then associated with their corresponding 100 × 100 km grid cell using the ‘join attributes by location’ tool. Environmental niche heterogeneity per grid cell was calculated as the mean of the standard deviations of all environmental variables extracted at species occurrences within a given grid cell. We tested for correlations between individual environmental variables and richness and speciation rate using multiple regressions in R (R Core Team, 2018). We initially conducted an ordinary least squares (OLS) regression and then to account for spatial autocorrelation we used a simultaneous autoregressive (SAR) lag model (Kissling & Carl, 2007) implemented using the ‘spdep’ package (Bivand & Piras, 2015). For each region, we used the number of nearest neighbours that optimized the fit of the model based on Akaike’s information criterion (AIC). We compared the SAR versus OLS model with AIC. In the linear regressions, we used the R2 (or Nagelkerke pseudo‐R2 in the case of the SAR; Nagelkerke, 1991) value from the best fitting model. Finally, we estimated the relative contribution of each environmental variable to global fern species richness using a multiple regression and accounting for spatial autocorrelation.
2 RESULTS
Species richness per grid cell ranged from 3 to 929 species (Figure 1). Unsurprisingly, we recovered a positive correlation between the number of collections in a grid cell and species richness (see Appendix S1.1). Our hotspot analysis identified eight distinct fern diversity hotspots based on species richness and endemism: Mesoamerica, the Greater Antilles, the tropical Andes, Southeastern Brazil, the Guianas, Madagascar, Malesia, East Asia (Figure 1; Figure S1.2). Within circumscribed hotspots, species richness ranged from 539 species in Madagascar to 2001 species in the tropical Andes; on a per‐area basis, richness was lowest in East Asia (7 spp. per 10,000 km2) and highest in the Greater Antilles (49 spp. per 10,000 km2) (Figure 2; Figure S1.3a). Within hotspots, richness per unit area was always higher at high elevations (> 1000 m) compared to low elevations (< 1000 m), ranging from 1.3× higher in the tropical Andes to 15.1× higher in the Greater Antilles (see Appendix S2). We recovered a mid‐elevation species richness peak on a global scale, with a slight shift to upper elevations in the tropics (Figure 3a). Furthermore, richness was generally highest at low latitudes (Figure 3b). Across hotspots, area was positively correlated with species richness (R2 = 0.93, p = 0.00001), however not significantly with the number of endemic species (R2 = 0.47, p = 0.06) (Figure S1.3a,b).



Endemism was higher in all hotspots compared to neighbouring non‐hotspots regions (Figure S1.2). A total of 6099 distinct species were recovered within all seven hotspots (Appendix S2). Percent endemism within hotspots ranged from 10% in the Guianas to 75% in Malesia (Figure 2). The proportion of endemic species was higher in high‐elevation subregions compared to low‐elevation subregions in all sites excluding Mesoamerica (high elevation: 27%, low elevation: 33%); the tropical Andes (high elevation: 38%, low elevation: 26%), Southeastern Brazil (high elevation: 32%, low elevation: 28%), Greater Antilles (high elevation: 30%, low elevation: 28%), Guianas (high elevation: 11%, low elevation: 10%), East Asia (high elevation: 72%, low elevation: 63%), Malesia (high elevation: 72%, low elevation: 72%) and Madagascar (high elevation: 77%, low elevation: 64%). Our CANAPE analyses indicated generally high levels of endemism within hotspots, but the type of PE differed between these regions (Figure S1.4). We found neo‐endemic grid cells in Madagascar, East Asia and sparsely in the tropical Andes (Figure S1.4). Conversely, palaeo‐endemic grid cells were mostly restricted to Malesia, with a smaller number occurring in Madagascar, as well as in non‐hotspot areas in Oceania and temperate South America (Figure S1.4). Areas of mixed (significant neo‐ and palaeo‐endemism) and super‐endemism (highly significant neo‐ and palaeo‐endemism) were widespread, appearing most conspicuously in Malesia, East Asia, Madagascar and the Greater Antilles, and some non‐hotspot areas in southern South America, southern Africa and northern Mexico (Figure S1.4).
Species present in the phylogeny were relatively evenly distributed across all fern families (Appendix S2). Tip speciation rates ranged from 0.005 to 0.631 events MY−1 and were highest in the families: Athyriaceae, Dryopteridaceae and Polypodiaceae. Average grid cell speciation rates ranged from 0.1 to 0.42 events MY−1, with speciation rates generally highest at low latitudes (Figures S1.5 and S1.6). These rates also varied between hotspots, with some of the highest rates found in the Guianas, tropical Andes and East Asia (Figure S1.7). Average grid cell speciation rates were also higher within hotspots compared to neighbouring non‐hotspot areas, excluding tropical Asia (Figure S1.8). Species assemblages in centres of neo‐endemism exhibited higher mean speciation rates than palaeo‐endemic centres (p < 0.001) (Figure 4a). Niche heterogeneity did not differ between centres of neo‐endemism and palaeo‐endemism, but neo‐endemism centres were more topographically heterogeneous compared to palaeo‐endemic centres (p < 0.001) (Figure 4b,c).

We found that hotspots harbour a disproportionate amount of temperature–precipitation climate space compared to the land area that they occupy (Figure 5). The hotspots with the largest climate space occupancy were the tropical Andes (52%), East Asia (33%), Malesia (27%) and Mesoamerica (26%) (Figure 5). The multiple regression of ecological variables by richness on a global scale demonstrated that variation in annual precipitation, coldest minimum temperature and soil carbon availability are most strongly predictive of global species richness patterns (Appendix S2). With varying correlation coefficients, species richness was positively correlated with niche heterogeneity in all sites, as well as all non‐hotspot sites (Figure S1.9). Species richness in hotspots was positively correlated with the percent of global climate space occupied (Figure S1.10). Environmental niche heterogeneity was higher in hotspots compared to adjacent non‐hotspot areas (Figure 6 inset), but the relationship between niche heterogeneity and speciation rate varied by region. In the Neotropics, we recovered a significant positive correlation between niche heterogeneity and speciation rates in all hotspots except Southeastern Brazil and the Greater Antilles (Figure 6). In the Palaeotropics, we found a significant, positive relationship in all hotspots (Figure 6).


3 DISCUSSION
The skewed distribution of life on Earth has attracted the attention of biologists for centuries (Humboldt & Bonpland, 1805). A growing body of literature increasingly supports the long‐held observation that mountains are home to some of the world’s most diverse ecosystems (Hoorn et al., 2013; Rahbek, Borregaard, Antonelli, et al., 2019); however, there is a lack of consensus on which factors drive this pattern (Perrigo et al., 2020; Rahbek, Borregaard, Colwell, et al., 2019). Here, we discuss our principal findings on global fern diversity, with the aims of understanding the evolution of this major plant lineage and providing insights on broader questions about the formation of global biodiversity hotspots.
3.1 Tropical mountains influence global patterns of species richness
Mountains, especially in the tropics, harbour disproportionately high levels of species richness and endemism relative to the land area they occupy (Rahbek, Borregaard, Colwell, et al., 2019). We demonstrate that ferns follow this pattern with striking prominence: 58% of all species occur in eight principally montane tropical and subtropical hotspots which together comprise just 7% of Earth’s land area (Figure 2; Appendix S2). This exceptional richness is mirrored by increased climate space in these regions. For example, the tropical Andes harbour more than 50% of global temperature–precipitation variation, despite covering just 1.2% of Earth’s land surface (Figure 5). Tropical and subtropical mountains with steep elevational gradients lead to greater ecological niche space compared to surrounding lowlands with lower fern diversity, and we recover a scalable, nearly universal, pattern of local fern species richness being positively correlated with available niche space (Figures S1.9 and S1.10). This is compounded by the fact that the volume of mountainous areas is much greater than the planimetric surface area they occupy. The pattern of outsized climate space occupancy is even more pronounced when only montane portions of these areas are considered: hotspot areas above 1000 m in elevation cover < 2.0% of Earth’s land surface but harbour over 45% of Earth’s fern diversity (Figure 2).
Combining all records (i.e., in hotspot and non‐hotspots), we find a strong latitudinal diversity gradient (Figure 3b) and a middle to upper elevational peak in species richness (Figure 3a), as reported in previous transect‐based studies (Hernández‐Rojas et al., 2020; Kessler et al., 2011; Kluge et al., 2006; Weigand et al., 2020). While the mid‐elevation richness peak is a well‐established phenomenon observed in many tropical plant and animal lineages (Rahbek, 1995), the pattern is particularly strong and consistent in ferns (Cardelús et al., 2006; Kessler et al., 2011; Watkins et al., 2006). Contention about the factors underlying the general mid‐elevation richness peak aside (Colwell & Lees, 2000), understanding why ferns adhere to this pattern so strongly is critical to explain variation in the group’s richness at broad scales.
We propose that the exceptional fern species richness in tropical middle to upper elevation forests (Figure 3a) is largely due to the prominence of many epiphytic fern lineages in these regions. Although epiphytic ferns do not appear to diversify more rapidly than terrestrial species (Schuettpelz & Pryer, 2009; Testo & Sundue, 2018), nearly perhumid conditions, due to significant water input from ‘horizontal precipitation’ (Vogelmann, 1973), allow for the persistence of a diverse vascular epiphyte flora in tropical cloud forests (Cardelús et al., 2006; Küper et al., 2004). Taxonomic representation among vascular epiphytes is highly skewed with a few successful lineages such as orchids and bromeliads comprising the majority of diversity (Gentry & Dodson, 1987). Along with orchids, ‘ferns are the predominant vascular epiphytes nearly everywhere’ (Gentry & Dodson, 1987), and can comprise up to 25% of the vascular epiphyte flora in some tropical cloud forests (Cardelús et al., 2006), despite making up less than 3% of global vascular plant diversity (PPG, 2016). In many middle to upper elevation tropical forests, epiphytic ferns are more diverse than their terrestrial counterparts; for instance in the cloud forests of Costa Rica epiphytic species can comprise 84% of local fern diversity (Watkins et al., 2006). The success of ferns in epiphytic niches is also clear when viewed from an evolutionary perspective: several of the most diverse fern clades (e.g. Asplenium, Elaphoglossum, Hymenophyllaceae, Polypodiaceae) are principally epiphytic (Schuettpelz & Pryer, 2009; Testo & Sundue, 2016). Though epiphytic ferns occur across a range of elevations, and even into temperate regions, their exceptional diversity in montane forests is clearly a principal component of the overall prominence of ferns in middle to upper elevation sites in the tropics.
3.2 Geographical patterns and drivers of speciation
On a global scale, we find that rapidly speciating lineages are primarily found in middle to lower latitudes (Figure S1.6), particularly in montane sites in East Asia and northern South America (Figure S1.5). Speciation rates of lineage assemblages in most regions are higher within hotspot grid cells compared to those in adjacent non‐hotspot grid cells (Figure S1.8). Two hotspots harbouring some of the most rapidly speciating lineages are the tropical Andes and East Asia (Figures S1.5 and S1.7). These sites correspond closely to mountain systems that have undergone recent and rapid orogeny and include some of the steepest elevational and climatic gradients on Earth (Clark et al., 2005; Hoorn et al., 2010). In addition to ferns (Figure S1.5), these same mountains also host other rapidly diversifying lineages (Madriñán et al., 2013; Xing & Ree, 2017). Recent studies have demonstrated a correlation between rates of mountain uplift and the exceptional speciation in plant groups such as bellflowers (Lagomarsino et al., 2016), bamboos (Ye et al., 2019), clubmosses (Testo et al., 2019) and orchids (Pérez‐Escobar et al., 2017). We find the tropical Andes and Hengduan Mountains include some of the youngest and most rapidly speciating fern lineages, such as Serpocaulon and Elaphoglossum in the Andes, and Athyrium and Lepisorus in the Hengduan Mountains (Appendix S2). However, we also find examples of rapid speciation outside of our circumscribed hotspots: both the Cape Region of Southern Africa and parts of Australia harbour rapidly diversifying fern assemblages (Figure S1.5). In these cases, it appears that, although there is evidence of recent radiations, overall fern diversity remains modest (Figure 1).
Within hotspots, we found a strong relationship between speciation rates of lineage assemblages and environmental niche heterogeneity (Figure 6); this is most evident in the tropical Andes, Mesoamerica and East Asia. Globally, we find sites of rapid speciation are generally in topographically complex regions (Figure S1.5). Geological, geographical and climatic factors underlie processes that play an important role in promoting speciation, including orogeny‐induced generation of novel niche space, rapid climatic turnover along steep elevational gradients, and local variation in soil and bedrock types (Condamine et al., 2018; Givnish et al., 2014; Hoorn et al., 2010, 2013, 2018; Rahbek, Borregaard, Antonelli, et al., 2019).
3.3 Factors influencing patterns of endemism
Both endemism and species richness are higher in montane sites when considered on a per‐unit‐area basis (Figure 2). However, total numbers of endemic species in most hotspots were similar in montane and lowland habitats (Figure 2); this contrasts sharply with a general pattern of increased endemism at high elevations that has been reported for many groups of organisms (Rahbek, Borregaard, Colwell, et al., 2019). Reasons for this departure remain uncertain, but we propose several possibilities. First, the high dispersibility of fern spores may allow species inhabiting patchily distributed montane habitats to maintain sufficient gene flow between populations to avoid isolation on ‘sky islands’ (Barrington, 1993), which is thought to be a principal mechanism by which narrowly endemic species form in other montane groups (Hughes & Eastwood, 2006). Second, available evidence suggests that fern lineages tend to show strong affinity to soil types, and species distributions at local scales may largely be mediated by edaphic features (Lehtonen et al., 2017; Tryon, 1972). Because diversity of soil and bedrock types is largely independent of elevation at broad spatial scales (Antonelli et al., 2018; Rahbek, Borregaard, Antonelli, et al., 2019), edaphic constraints on species distributions may help explain the patterns observed here. Historical climate factors not evaluated in our analyses may have played an important role in shaping patterns of endemism along elevational gradients. As hypothesized by Kluge et al. (2006) for Mesoamerican ferns, upwards shifts of forest assemblages during periods of global warming in the Quaternary may have reduced the number of high‐elevation endemics through local extinctions on isolated mountaintops.
Although levels of endemism are high across global hotspots (Figure S1.2), the type of PE varies by site (Figure S1.4). Parts of the tropical Andes, Madagascar and East Asia include sites with high levels of neo‐endemism (Figure S1.4); Malesia, in contrast, harbours most of the world’s centres of palaeo‐endemism and mixed‐endemism (areas with both high neo‐endemism and palaeo‐endemism; Figure S1.4). The distribution of these endemism centres provides insight into the disparate evolutionary histories of ferns in these regions. Within the tropical Andes and the Hengduan Mountains, evidence suggests that areas with incidence of neo‐endemism are due to recent in situ speciation at high elevations (Figures S1.5 and S1.7; Sánchez‐Baracaldo & Thomas, 2014; Wei et al., 2018; Xing & Ree, 2017). As a rule, these habitats are geologically complex and distribution of their vegetation is strongly impacted by climatic variation along an elevation cline (Clark et al., 2005; Hoorn et al., 2010). Dynamics between climatic variation across, and seasonal stability within, elevational transects in these tropical and subtropical mountains is hypothesized to reduce the range size of many groups of organisms leading to high allopatric speciation (Janzen, 1967; Polato et al., 2018). Further study is needed to better understand how ‘sky islands’ and climate‐driven shifts in their extent impacts patterns of speciation and endemism in ferns and other groups; methods that incorporate spatially explicit models of habitat connectivity through time (Flantua et al., 2019) may prove to be particularly insightful.
In contrast, the predominance of palaeo‐endemism in Malesia highlights the importance of isolation and long‐term climatic stability in maintaining the diverse and unique species assemblages found in this hotspot (Figure S1.4). Although we find evidence of rapid, recent speciation in some areas (e.g. the Papuan Highlands and Borneo; Figure S1.5), Malesia has the lowest mean speciation rate of any of our hotspots (Figure S1.7). Compared to other global hotspots, the high fern diversity observed in Malesia appears to be largely due to the persistence of ancient endemic lineages. Most Malesian ferns are found nowhere else on Earth (Figure 2), including members of relictual, once widespread, families such as Cystodiaceae and Matoniaceae, which are represented in the fossil record dating back hundreds of millions of years (Choo & Escapa, 2018; Regalado et al., 2018). The importance of Malesia as both a biodiversity hotspot and ‘museum’ has been noted for other groups (Steenis, 1979) and reflects the signature of the region’s complex palaeogeographical and climatic histories on its biota. Since the Cretaceous, Malesia has experienced a complex history of land submersion and emergence events, resulting in a continually changing mosaic island chains with varying levels of connectivity (Lomolino et al., 2006). This highly dynamic landscape history led to a multitude of local speciation, extinction and migration events that resulted in conspicuous patterns of distributions observed today, and have contributed to exceptional diversification bursts in numerous lineages (Moyle et al., 2016; Thomas et al., 2012). At the same time, Malesia has acted as a refuge for many relictual lineages that were formerly widespread, suggesting that this region has experienced reduced rates of extinction compared to continental hotspots. The persistence of ancient groups on these and other islands (such as Madagascar) is likely aided by the thermal stability of insular systems and lower rates of biotic turnover (Cronk, 1997; Fjeldså et al., 2012).
3.4 Conclusions, caveats and future directions
Recent work focusing on tropical mountains has highlighted the geologic and climatic factors that help generate and maintain their exceptional biodiversity, but most broad‐scale analyses have focused on a few well‐known groups, such as mammals, birds and amphibians (Antonelli et al., 2018; Perrigo et al., 2020; Quintero & Jetz, 2018; Rahbek, Borregaard, Antonelli, et al., 2019). In this study, we close a significant gap in our knowledge with a global analysis of patterns of diversity in a major plant lineage: ferns. By integrating phylogenetic, ecological and georeferenced data, we demonstrate the following points. By our estimates, 58% of global fern diversity is concentrated in only eight biodiversity hotspots and nearly half (~45%) is concentrated in montane sites (Figures 1 and 2). Second, these hotspots occupy a disproportionate amount of climatic and environmental niche space compared to their total land area, which acts as a strong predictor of its species richness (Figure 5; Figures S1.9 and S1.10). Third, species assemblages with the highest rates of speciation are generally found at middle to lower latitudes and higher elevations (Figure S1.6), especially in areas with high environmental heterogeneity (Figure 6). Finally, patterns of endemism vary across hotspots: neo‐endemism is prevalent in areas with high, young mountains, whereas tropical islands with stable climates and relatively low biotic turnover, act as refugia for ancient lineages (Figure 3). Taken together, these findings suggest that mountain building and associated formation of heterogeneous environments promote rapid species proliferation through reduced dispersal and gene flow, while long‐term ecological and climatic stability helps maintain global diversity in this ancient clade by buffering against extinction in formerly widespread groups.
There are several limitations to the present study. First, it is well known that quality of available species occurrence data is highly variable and that misidentified records or those with problematic localities can bias biodiversity analyses (Zizka et al., 2019). Although we made extensive efforts to ensure the quality of the species occurrence data used through a series of automated and manual cleaning steps—resulting in the removal of 50% of the records in the original dataset—some problematic records almost certainly remain in our final dataset. In addition, political and geographical biases in data availability are evident. For example, few records of ferns from India are available on GBIF, despite the fact that its fern flora is diverse and relatively well‐studied (Fraser‐Jenkins et al., 2017). Similarly, the marked difference in record density and estimated species richness between the Indonesian and Papuan sides of New Guinea appear to reflect differences in collection effort and digitization, rather than true diversity patterns (Cámara‐Leret et al., 2020). Unsurprisingly, collection effort and inferred species richness are positively correlated globally (Figure S1.1); nonetheless, our estimates of species richness are similar to species counts obtained from floristic treatments for most regions (Appendix S2). These challenges are shared by other studies relying on global species occurrence datasets (Yesson et al., 2007) and interpretation of our findings should be made with appropriate caution. Our findings are also shaped in part by the phylogeny used in several analyses, ranging from speciation rate estimation to characterization of phylogenetic endemism centres. While this phylogeny enabled us to maximize the number of taxa with evolutionary data, more than half of known fern species are yet to be placed in a phylogeny, and some genera remain unsampled altogether (Appendix S2). In addition, although the phylogeny used in this study is largely topologically concordant with those of other studies (Lehtonen et al., 2017; PPG, 2016; Schuettpelz & Pryer, 2009), uncertainty in relationships among lineages and with respect to divergence times could impact some of the findings reported here. Together with improved species occurrence data, phylogenies with improved taxon sampling, increased numbers of loci and more sophisticated divergence time estimation methodologies should provide improved resolution of the patterns explored here.
With these caveats considered, we suspect that our main conclusions are valid. The hotspots of richness and endemism that we recover correspond closely to those reported previously in floristic studies as well as plot‐based analyses examining fern diversity at local to global scales (Hernández‐Rojas et al., 2020; Kessler et al., 2011; Kluge et al., 2006; Moran, 1995, 2008; Nagalingum et al., 2015; Tryon, 1972; Weigand et al., 2020; Xing & Ree, 2017). In addition, our findings are similar to those found for other taxonomic groups and with varying methodological approaches (Antonelli et al., 2018; Myers et al., 2000; Perrigo et al., 2020; Rahbek, Borregaard, Antonelli, et al., 2019). In light of this, while we anticipate that insight into global fern diversity will improve as more data become available for currently underrepresented regions, we see little reason to expect significant departures from the broad‐scale patterns reported here.
We anticipate that future work will continue to investigate drivers of species diversity by leveraging methodological advances and more robust datasets. In particular, generating a comprehensive understanding of the importance of mountains for biodiversity requires improved integration of biology and geology. Though our findings add to a substantial body of work highlighting the exceptional richness, endemism and rapid diversification of montane groups of organisms, we did not attempt to assess how the timing and rates of mountain uplift events correspond with patterns of lineage diversification. Environment‐dependent models of diversification allow for such correlations to be tested (Condamine et al., 2018), but robust estimates of uplift rates are lacking for many mountain systems and applying these models to a global analysis remains challenging. Similar models can be extended to better integrate species diversification estimates with other important historical factors, including changes in global climate and shifts in the presence and connectivity of islands, both oceanic islands and high‐elevation ‘sky islands’. Integrating these models with improved species‐level data will prove particularly valuable as we seek to better understand what factors influence species distributions, patterns of niche evolution and lineages’ evolutionary trajectories through time.
Acknowledgements
We thank David Barrington and Michael Kessler for conversations on plant biogeography and Morgan Southgate for assistance on spatial analyses. We are indebted to Michael Hassler for sharing data from the Checklist of Ferns and Lycophytes of the World, and Shawn Laffan for support with the Biodiverse analysis package. Finally, we wish to recognize that this work would not be possible without the efforts of countless collectors, taxonomists, curators and herbarium staff.
References
Biosketches
Jacob Suissa is a PhD candidate in the Department of Organismic and Evolutionary Biology at Harvard University, as well as a Fellow of the Arnold Arboretum. His research interests include the evolution of ferns and lycophytes at both the organismal and macroevolutionary scale.
Michael Sundue is a research associate professor and curator of the Pringle Herbarium at the University of Vermont, USA. His research interests include palaeobotany and macroevolution.
Weston Testo is a postdoctoral research associate at the University of Gothenburg, Sweden. He is interested in the evolution of fern and lycophyte diversity, especially in the American tropics.
Author contributions: JSS and MAS initially conceived the project. JSS, MAS and WLT developed the project to its final state. JSS and WLT conducted the research and created figures. JSS, WLT and MAS wrote the manuscript.